单细胞转录组基础分析三：降维与聚类

单细胞测序技术是近年最大的生命科学突破之一，相关文章频繁发表于各大顶级期刊，然而单细胞数据的分析依然是大家普遍面临的障碍。本专题将针对10X Genomics单细胞转录组数据演示各种主流分析，包括基于Seurat的基础分析、以及基于clusterProfiler、Monocle、SingleR等R包的延伸分析。不足之处请大家批评指正，欢迎添加Kinesin微信交流探讨！（扫描文末二维码）

单细胞研究的一个重要内容就是对细胞进行分类和识别，细胞分类的基础是它们的基因表达特征。然而人类有3万多个基因，1个基因的值代表了定义细胞表达特征的1个维度，3万多个基因就是3万多个维度的超高维度空间，这是人类无法理解的数据。Seurat会用三个步骤把细胞的表达特征在我们习惯的二维空间呈现出来：首先是选择2000个（默认）表达值变化大的基因代表细胞转录谱，其次使用主成分分析将2000维的信息投射到50个维度，并提取前10-20个维度（人工选择）的信息代表细胞的转录特征，最后使用非线性降维方法（tSNE或UMAP）将这10-20个PC值降维到二维空间。经过上述操作，转录谱的特征信息会损失一些，但是大部分转录特征会在二维空间呈现出来。借用一个使用PCA算法压缩图片的例子来展示降维操作的效果：

第一张图是有37270个颜色维度的原图，后面是提取不同数量的主成分的重构图，可以发现少量的主成分就可以呈现原图的特征。

寻找高变基因

Seurat负责筛选高变基因的函数是FindVariableFeatures()，它并不删除scRNA对象中的非高变基因。找出的结果可以通过VariableFeatures()函数获取。

library(Seurat)
library(tidyverse)
library(patchwork)
rm(list=ls())
dir.create("cluster")
scRNA <- readRDS("scRNA.rds")
scRNA <- FindVariableFeatures(scRNA, selection.method = "vst", nfeatures = 2000) 
top10 <- head(VariableFeatures(scRNA), 10) 
plot1 <- VariableFeaturePlot(scRNA) 
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) 
plot <- CombinePlots(plots = list(plot1, plot2),legend="bottom") 
ggsave("cluster/VariableFeatures.pdf", plot = plot, width = 8, height = 6) 
ggsave("cluster/VariableFeatures.png", plot = plot, width = 8, height = 6)

选择的高变基因可通过cluster/VariableFeatures.png文件查看，横坐标是某基因在所有细胞中的平均表达值，纵坐标是此基因的方差。红色的点是被选中的高变基因，黑色的点是未被选中的基因，变异程度最高的10个基因在如图中标注了基因名称。

在进行PCA降维之前还有要对数据进行中心化（必选）和细胞周期回归分析（可选），它们可以使用ScaleData()函数一步完成。这两项分析虽然可以使用一个函数完成，但是我觉得有必要分开讲一下它们的原理与作用。

数据中心化

数据中心化是PCA分析和一些可视化步骤的必须要求，它使数据具有统一的尺度（scale），其运行效果示意图如下所示：

左边的图是原始数据画的散点图，右边的图是用scale之后的数据画的。细心的同学可能发现了，这两张图点的相对位置没变，但是坐标轴的尺度改变了。

Scaling代码如下：

##如果内存足够最好对所有基因进行中心化
scale.genes <-  rownames(scRNA)
scRNA <- ScaleData(scRNA, features = scale.genes)
##如果内存不够，可以只对高变基因进行标准化
scale.genes <-  VariableFeatures(scRNA)
scRNA <- ScaleData(scRNA, features = scale.genes)

scRNA对象中原始表达矩阵经过标准化和中心化之后，已经产生了三套基因表达数据，可以通过以下命令获得：

#原始表达矩阵
GetAssayData(scRNA,slot="counts",assay="RNA") 
#标准化之后的表达矩阵              
GetAssayData(scRNA,slot="data",assay="RNA")
#中心化之后的表达矩阵 
GetAssayData(scRNA,slot="scale.data",assay="RNA") 

细胞周期回归

上一步找到的高变基因，常常会包含一些细胞周期相关基因。它们会导致细胞聚类发生一定的偏移，即相同类型的细胞在聚类时会因为细胞周期的不同而分开。如在Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing 中 vCAFs和cCAFs中只有cell cycle genes的表达差异。

查看我们选择的高变基因中有哪些细胞周期相关基因：

CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(scRNA))

细胞周期评分

g2m_genes = cc.genes$g2m.genes
g2m_genes = CaseMatch(search = g2m_genes, match = rownames(scRNA))
s_genes = cc.genes$s.genes
s_genes = CaseMatch(search = s_genes, match = rownames(scRNA))
scRNA <- CellCycleScoring(object=scRNA,  g2m.features=g2m_genes,  s.features=s_genes)

以上代码运行之后会在[e-mail]中添加S.Score、G2M.Score和Phase三列有关细胞周期的信息。

查看细胞周期基因对细胞聚类的影响

scRNAa <- RunPCA(scRNA, features = c(s_genes, g2m_genes))
p <- DimPlot(scRNAa, reduction = "pca", group.by = "Phase")
ggsave("cluster/cellcycle_pca.png", p, width = 8, height = 6)

右图是seurat细胞周期回归分析教程的案例，提示这种情况需要剔除细胞周期对聚类的影响。左图是我们这次分析的pbmc的结果，看情况没必要做细胞周期回归分析。

##如果需要消除细胞周期的影响
#scRNAb <- ScaleData(scRNA, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(scRNA))

PCA降维并提取主成分

PCA降维

scRNA <- RunPCA(scRNA, features = VariableFeatures(scRNA)) 
plot1 <- DimPlot(scRNA, reduction = "pca", group.by="orig.ident") 
plot2 <- ElbowPlot(scRNA, ndims=20, reduction="pca") 
plotc <- plot1+plot2
ggsave("cluster/pca.pdf", plot = plotc, width = 8, height = 4) 
ggsave("cluster/pca.png", plot = plotc, width = 8, height = 4)

运行之后会在cluster文件夹下面得到pca图，如下所示：

左图是根据主成分1和2的值将细胞在平面上展示出来，右图展示前20个主成分的解释量（pca中等同于标准差）。后续分析要根据右图选择提取的pc轴数量，一般选择斜率平滑的点之前的所有pc轴，此图我的建议是选择前18个pc轴。

pc.num=1:18

获取PCA结果 

此部分代码为分析非必须代码，不建议运行！！！
#获取基因在pc轴上的投射值
Loadings(object = scRNA[["pca"]])
#获取各个细胞的pc值
Embeddings(object = scRNA[["pca"]])
#获取各pc轴解释量方差
Stdev(scRNA)
#查看决定pc值的top10基因， 此例查看pc1-pc5轴
print(scRNA[["pca"]], dims = 1:5, nfeatures = 10) 
#查看决定pc值的top10基因在500个细胞中的热图，此例查看pc1-pc9轴
DimHeatmap(scRNA, dims = 1:9, nfeatures=10, cells = 500, balanced = TRUE)

细胞聚类

此步利用 细胞-PC值 矩阵计算细胞之间的距离，然后利用距离矩阵来聚类。其中有两个参数需要人工选择，第一个是FindNeighbors()函数中的dims参数，需要指定哪些pc轴用于分析，选择依据是之前介绍的cluster/pca.png文件中的右图。第二个是FindClusters()函数中的resolution参数，需要指定0.1-0.9之间的一个数值，用于决定clusters的相对数量，数值越大cluters越多。

scRNA <- FindNeighbors(scRNA, dims = pc.num) 
scRNA <- FindClusters(scRNA, resolution = 0.5)
table([e-mail]$seurat_clusters)
metadata <- [e-mail]
cell_cluster <- data.frame(cell_ID=rownames(metadata), cluster_ID=metadata$seurat_clusters)
write.csv(cell_cluster,'cluster/cell_cluster.csv',row.names = F)

运行之后会在rstudio中显示每个cluster的细胞数量，cluter文件夹中也会得到细胞-cluter的映射文件cell_cluster.csv

非线性降维

此步也是利用 细胞-PC值 矩阵作为输入数据

#tSNE
scRNA = RunTSNE(scRNA, dims = pc.num)
embed_tsne <- Embeddings(scRNA, 'tsne')
write.csv(embed_tsne,'cluster/embed_tsne.csv')
plot1 = DimPlot(scRNA, reduction = "tsne") 
ggsave("cluster/tSNE.pdf", plot = plot1, width = 8, height = 7)
ggsave("cluster/tSNE.png", plot = plot1, width = 8, height = 7)
#UMAP
scRNA <- RunUMAP(scRNA, dims = pc.num)
embed_umap <- Embeddings(scRNA, 'umap')
write.csv(embed_umap,'cluster/embed_umap.csv') 
plot2 = DimPlot(scRNA, reduction = "umap") 
ggsave("cluster/UMAP.pdf", plot = plot2, width = 8, height = 7)
ggsave("cluster/UMAP.png", plot = plot2, width = 8, height = 7)
#合并tSNE与UMAP
plotc <- plot1+plot2+ plot_layout(guides = 'collect')
ggsave("cluster/tSNE_UMAP.pdf", plot = plotc, width = 10, height = 5)
ggsave("cluster/tSNE_UMAP.png", plot = plotc, width = 10, height = 5)
##保存数据
saveRDS(scRNA, file="scRNA.rds")

运行之后会得到TSNE图和UMAP图，以及细胞在这两张图上的坐标值文件embed_tsne.csv和embed_umap.csv

获取帮助

本教程的目的在于把常用的单细胞分析流程串起来，适合有一定R语言基础的朋友参考。分析原理和代码我没有详细解释，官网有详细的教程和权威的解释，翻译好的中文教程也有多个版本，有兴趣的可以搜索一下。如果您学习本教程有一定困难，可以分享此篇文章到朋友圈，截图微信发给Kinesin（文末二维码），我会抽时间组群直播讲解单细胞数据分析的全过程。本专题所用的软件、数据、代码脚本和中间结果，我打包放在了百度云上，需要的朋友添加Kinesin微信索取。